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^\ ' Abstract. We present a new approach to the statistical study and modelling of number counts of faint point sources in as- 

tronomical images, i.e. counts of sources whose flux falls below the detection limit of a survey. The approach is based on the 
' theory of a-stable distributions. We show that the non-Gaussian distribution of the intensity fluctuations produced by a generic 

^ I point source population - whose number counts follow a simple power law - belongs to the a-stable family of distributions. 

Even if source counts do not follow a simple power law, we show that the a-stable model is still useful in many astrophysical 
scenarios. With the ff-stable model it is possible to totally describe the non-Gaussian distribution with a few parameters which 
are closely related to the parameters describing the source counts, instead of an infinite number of moments. Using statisti- 
^— ^ ■ cal tools available in the signal processing literature, we show how to estimate these parameters in an easy and fast way. We 

' demonstrate that the model proves valid when applied to realistic point source number counts at microwave frequencies. In the 

, case of point extragalactic sources observed at CMB frecuencies, our technique is able to successfully fit the P(D) distribution 

of deflections and to precisely determine the main parameters which describe the number counts. In the case of the Planck 
mission, the relative errors on these parameters are small either at low and at high frequencies. We provide a way to deal with 
the presence of Gaussian noise in the data using the empirical characteristic function of the P(D). The formalism and methods 
Q ' here presented can be very useful also for experiments in other frequency ranges, e.g. X-ray or radio Astronomy. 
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1. Introduction galaxies whose typical angular size is much smaller than the 

observing beam width (hence the name 'point sources') are es- 
C3 The study of intensity (or temperature) fluctuations in the pecially diflicult to deal with. The galaxies that contribute to the 
Cosmic Microwave Background (CMB) radiation has be- observed total signal are very diverse, corresponding to objects 
come one of the milestones of modern cosmology due to at different redshifts and with different physical properties. This 
two main reasons. On the one hand, the angular power spec- makes it impossible to estabUsh a single spectral behaviour 
trum of the fluctuations allow us to place tight constraints on for all of them, thus hampering the performance of methods 
the fundamental cosmological parameters (see, e.g., Bennett that use multi-frequency observations to achieve an efficient 
et al. |2M3a). For a recent review on the study of CMB component separation. The spatial distribution of faint EPS is 
anisotropies, see Hu & Dodelson ( 2002 ). On the other hand, roughly uniform across flie sky even in presence of source clus- 
the study of the different physical sources (foregrounds) that tering. Therefore, Galactic cuts useful to avoid Galactic fore- 
contribute to the incoming radiation at microwave wavelengths grounds such as synchrotron, dust and free-free emission, are 
has a great scientific relevance in itself (De Zotti et al. 1999 1, not applicable to avoid EPS contamination. 
Therefore, a great deal of effort has been devoted to the task gPS contamination occurs when a set of point-like sources 
of separating the different components that are present in CMB ^jj^ intensities distributed foUowing a general law, usually 
maps. In general, component separation techniques take advan- modeUed as a power law, is observed by a detector using a 
tage of the statistical behaviours of each component to distin- gjygjj instrumental response. The final signal is a mixing in 
guish among them. Hence, it is important to have good statisti- ^^ich the brightest sources are still individually detectable over 
cal models for each of the components under study. ^ 'confusion noise' generated by the contributions of faint, un- 
Among the different foregrounds that appear in CMB ob- resolved sources; this situation is very common in astronomi- 
servations, extragalactic point sources (EPS), i.e. individual cal images and it has been studied first at radio and X-ray fre- 
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quencies. The intensity (or temperature) distribution given by 
unresolved sources is strongly non-Gaussian and shows long 
positive tails. This kind of behaviour is known in the signal 
processing literature as 'impulsive noise' . 

The effect of the confusion noise is two fold: on one 
hand the mean value of this noise is positive, producing a 
'source monopole' (integrated extragalactic background) that 
has to be summed up to the other components and, on the 
other hand, it gives rise to small scale intensity fluctuations. 
At microwave frequencies, the fluctuations generated by unde- 
tected sources can severely hamper the detection of true CMB 
anisotropics (Franceschini et al. 1989 1. Recently, Toffolatti et 
al. C1998 hereafter T98) presented a detailed analysis of the 
effect of point sources on CMB anisotropy maps. By exploit- 
ing a cosmological evolution model for radio and far-IR se- 
lected sources, they made precise predictions on source counts, 
on confusion noise and on the angular power spectrum due to 
undetected sources. In particular, they showed that the contri- 
bution of EPS will be very relevant at the lower and higher fre- 
quency channels of the future ESA Planck mission (Mandolesi 
et al. 119981 Puget et al. 1998 1. As for radio source counts, the 
predictions of T98 have been confirmed by the first year data of 
the NASA Wilkinson Microwave Anisotropy Probe (WMAP) 
mission (Bennett et al. 2003bl, al least up to frequencies of 
~ 30-40 GHz. The WMAP all sky cat alogue of bright ex- 
tragalactic sources (Bennett et al. I2003b> consists of 208 ob- 
jects with fluxes S > 0.9 H- 1 Jy on a sky area of 10.38 sr 
(\b"\ > 10°) whereas flie T98 model predicts 270-280 sources 
at 30 GHz in the same area with an average off^set of ~ 0.75 be- 
tween observed and predicted number of extragalactic sources. 
Moreover, the distribution of energy spectral indexes (i.e. a, 
S oc v^") of sources in the WMAP sample peaks around 
a - 0.0, which is exactly the mean value of the energy spec- 
tral index adopted in T98 for "flat"-spectrum sources, i.e. the 
dominant source population at these frequencies (the fraction 
of "steep'-spectrum sources being ~ 10 15%). It is also no- 
ticeable that the brightest source detected by WMAP has a flux 
density of 5 ^25 Jy which, again, corresponds to the flux value 
for which the T98 model predicts 1 source all over the sky. 
Another important result confirming the estimates of the T98 
cosmological evolution model is that a good agreement is cur- 
rently found between predictions based on that model and data 
on the excess angular power spectrum at small angular scales 
as well as on the angular bispectrum detected in the WMAP Q 
and V bands (Komatsu et al. 120031 Argiieso et al. "2003). Two 
other independent samples of extragalactic sources at 31 and 
34 GHz - from CBI (Mason et al. l2003t and VSA (Taylor et 
al. '2003 ) experiments, respectively - show that the T98 model 
correctly predicts number counts down to, at least, 5^10 mJy. 
Therefore, we can confidently use the T98 model for simulating 
Poisson distributed EPS in CMB sky maps, at least up to 50- 
100 GHz. At higher frequencies, more recent models can give 
a better fit to current data on source counts (see section 15. 1.2> . 
As the outcomes of the methods to be presented here are model 
independent, we still use the T98 model throughout the paper 

The current low sensitivity of detectors at CMB frequen- 
cies makes it impossible to test directly model counts down to 
fainter fluxes. On the other hand, more information on counts 



of faint sources, i.e. sources with fluxes fainter than the de- 
tection threshold of a given experiment, can be extracted by 
the analysis of the intensity fluctuations of point sources. The 
probability density function pdf of fluctuations due to unde- 
tected point sources, as a first step to the modelling of the 
confusion noise, has been studied since the middle of the last 
century (Scheuer [T957l [1974. Condon [T974I Hewish [7961^ . 
These works have shown that it is possible to find analytical ex- 
pressions for the characteristic function of the pdf (in Fourier 
space), but not for the pdf in the real space. This fact has ham- 
pered the development of specific statistical tools to deal with 
the EPS confusion noise. 

In analysing a sky map there are two traditional ways of de- 
termining the main statistical properties of a given source popu- 
lation. One possibility is to detect the brightest point sources in 
a given data set, e.g. using a linear filter to detect them, and then 
obtain parameters such as the number counts, their slope, etc. 
For example, Vielva et al. J2003t detect point sources in realis- 
tic Planck simulations using a Mexican Hat Wavelet technique 
and compare the number of detections with the input number 
counts, which correspond to the T98 model. The other possi- 
bility is to directly study the pdf of the confusion noise which, 
in general, is mixed with the signal coming from CMB and 
the other foregrounds plus instrumental noise. This is generally 
performed using statistical indicators such as the moments up 
to a certain degree (see, e.g., Rubino-Martm & Sunyaev 2003] 
and Pierpaoli 2003 1 or the non-Gaussianity of the wings. A 
computationally more complex way is to calculate numerically 
the theoretical pdf assuming some model for source counts and 
trying to fit it to the data (Condon & Dressel 1978 Franceschini 
et al. 1989^ Anyway, the lack of an analytical form for the pdf 
makes difficult to establish the optimal estimator of its param- 
eters. In particular, it is not clear how many moments are nec- 
essary to characterise the pdf (in principle, infinite of them) or 
which ones are more appropriate for extracting information (it 
is generally assumed, on the basis of mere intuition, that the 
third order moment should be one of them). 

In this paper we will focus on the application of a novel for- 
malism, the a-stable distributions, to model the pdf of the in- 
tensity fluctuations due to point extragalactic sources, a-stable 
distributions are known to be very efficient in modelling im- 
pulsive noise. They have a number of interesting mathemati- 
cal properties that make them very attractive; in particular, it 
is possible to show that the Gaussian distribution is a special 
case of the more general class of a-stable distributions and that 
ff-stable distributions satisfy a generalised form of the central 
limit theorem. Moreover, in this work we show that the pdf of 
a theoretical power law representing the number counts of ex- 
tragalactic sources observed with a fiUed-aperture instrument 
must follow exactly an ff-stable distribution.' The great advan- 



' The same analysis can also be applied to images obtained by in- 
terferometry techniques (insensitive to the zero level). In this case, the 
resulting pdf shows the same shape as for filled-aperture images but 
it is centred at 0. If filled-aperture measurements are currently per- 
formed by means of dual-beam scans, subtracting one signal form the 
other, filled-aperture and interferometry techniques yield similar pdfs 
(see, e.g., De Zotti et al. .l996.) . 
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tage is that o'-stable distributions are completely described by a 
small number of parameters, what makes them relatively easy 
to deal with. Optimal techniques already existent in the sig- 
nal processing literature are easy to adapt to directly extract 
the main parameters of the source number counts (namely, the 
slope of the number counts power law and its normalisation) 
without having to resort to clumsy statistics. Finally, the meth- 
ods can be generalised for dealing with mixtures of signals, as 
is the case when the EPS population is added to Gaussian in- 
strumental noise. 

Therefore, this work aims at four objectives: first, we will 
demonstrate that a generic extragalactic point source popula- 
tion whose number counts follow a power law inevitably leads, 
when observed with a filled-aperture instrument, to a pdf which 
belongs to the family of c-stable distributions. Then, we will 
devote some time to introduce the c-stable distributions to the 
community, reviewing their main properties and giving the nec- 
essary references for further reading. In a third part of this pa- 
per we will discuss how a-stable distributions can be used to 
model and study more realistic cases, such as the truncated 
power law, and we will show that the model is a good ap- 
proximation for most cases of astronomical interest. Finally, 
we will discuss the case in which other astronomical or instru- 
mental signals -such as the Cosmic Microwave Background 
(CMB) radiation or instrumental noise- are mixed with the 
point sources. In that case, a method will be suggested that is 
able to obtain the parameters of the point source deflection dis- 
tribution even in presence of 'contamination'. 

The structure of this paper is as follows: in Sect. |2lwe re- 
view the basics of the derivation of the characteristic function 
of the deflection distribution. In Sect. |3] the a-stable distribu- 
tions and their main properties are introduced. Section0]deals 
with the extraction of the physical parameters of the EPS pop- 
ulation using the a-stable formalism. In Sect. |5] we study the 
application of the formalism to more realistic source models, 
using as an example the T98 point source model. Parameter 
estimation of (T-stable processes mixed with Gaussian noise is 
considered in section|6l A few considerations about the imple- 
mentation of these techniques for the future Planck mission are 
given in sectional Finally, in Sect. |8] we summarise our con- 
clusions. 



2. Source counts and the deflection probability 
function P(D) 

Let us consider a population of EPS whose difl^erential number 
counts can be described in a power law form: 



n(S)^kS-\ S >0, 



(1) 



where t] is the slope of the diff'erential counts power law, k is 
called its normalisation and S is the intrinsic flux. The sources 
are assumed to be distributed uniformly across the sky and, at 
the moment, we will assume that eq. Q holds for all S > 0. 
The sources are now observed with an instrument whose angu- 
lar response is f{6, (p), not necessarily normahsed to unity at the 



peak. Then, the mean number of sources responses of intensity 
X = f(6, (p)S in the beam at any time is 



R(x) 



m<p) 



dO. 



Substituting eq. Q into eq. we have that 

Rix) = kQex''', 

where 

a = J [mcf>)]''-'dQ. 



(2) 



(3) 



(4) 



is a geometrical factor called effective beam solid angle. Let 
us now define the deflection D as the fluctuation field that is 
observed, that is D = / - (/), where / is the intensity at a given 
point (time) and (/) is its average value, i.e. (/) represents the 
extragalactic background due to undetected EPS. Let us define 
the characteristic function of a given function g{x) as 



G(w) 



f 

•J —cx 



8(x)' 



"dx. 



(5) 



Scheuer ( 1957 1 showed that the characteristic function of the 
probability distribution P{D) is related to the characteristic 
function of R(x) through 



iff(w) = exp [r(w) - r(Q)] , 



(6) 



where ilr(w) and r(w) are the characteristic functions of P(D) 
and R{x), respectively. 

Using eqs. Q, and (|6j it is possible to calculate the char- 
acteristic function i/'(w). This calculation has been performed 



by several authors, including Scheuer J1957> . Condon J 19741 



Barcons ( 1992) and Franceschini et al. ( I1989t . After some ef- 
fort we obtain 



i^(w) - exp < ifiw - y |w|' 



iaTi\ 

1 + iy?sgn(w)tanl — 1 



(7) 



where the parameters o-, y and p relate to the physical pa- 
rameters of the EPS and of the detector through 



1, 



^, 1 -I- orV / 1 - a 



r(^)r(^)sin(f)' 



lim a 



\ - a £nO+ 



(8) 
(9) 

(10) 
(11) 



The terms in eq. have been arranged in this way for a 
reason that will be clear in the following. The second equal- 
ity in eq. (|5Jl is due to the properties of the gamma func- 
tion (Abramowitz "1971 1, and as far has we know it hasn't 
been noticed before now. The previous equations are valid for 
1 < 77 < 3. For 77 > 2 the parameter yu is not finite, a situation 
equivalent to the classic Olbers' paradox in which the observed 
integrated flux density is infinite in aU directions of the sky. 
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Equation Q has an important drawback: to obtain the pdf 
of the deflections, PiD), it is necessary to make the inverse 
Fourier transform of iff{w) which, in general, cannot be eval- 
uated analytically. Although it can be performed numerically, 
the computational cost can be high if many different realisa- 
tions are needed for a particular task, and numerical integration 
does not lead to closed form solutions. Instead of doing that, let 
us see what can be learnt from the characteristic function itself. 
As we will see in the next section, eq. corresponds exactly 
to the actual definition of a family of distributions that is known 
in the statistical signal processing literature as a-stable distri- 
butions. 



3. A short introduction to a-stable distributions 

In this section we will make a summary introduction to the a- 
stable distributions, in order to clarify some concepts that will 
be used later. Starting from the particular case of the EPS P{D) 
characteristic function in eq. Q, that we will see that belongs 
to the a-stable family, we will proceed to a more general de- 
scription of a-stable phenomena that could be useful for as- 
tronomers. 

The characteristic function in eq. corresponds to a pdf 
that in general should be calculated numerically and that ex- 
hibits heavier tails than a Gaussian distribution. In Fig. the 
pdf corresponding to eq. (0 with [i - 1, y = 1 and three dif- 
ferent values of a (1.1, 1.5 and 1.9) are shown. The presence 
of heavy tails mean that 'glitches' are more likely to occur than 
in the Gaussian case. In the signal processing literature, proba- 
bility density functions with tails heavier than the Gaussian are 
called to be impulsive. A process is impulsive if it takes large 
values that significantly deviates from the mean value with non- 
negligible probability. These large values often appear as con- 
spicuous outlayers. Impulsive processes are ubiquitous in many 
'real-world' problems, from atmospheric noise caused by elec- 
tric discharges to financial time series data. For more informa- 
tion on impulsive noise, see Kuruoglu ( 1998 1. A great deal of 
effort has been done to model impulsive processes; is in that 
context that the a-stable distributions have experienced popu- 
larity. Other models that deal with impulsive processes, such 
as the Middleton's (Middleton 1977 1, Cauchy and Student-T 
models, are very case specific while the a-stable model is gen- 
eral and has a strong theoretical justification. 

Although a-stable distributions were known from the be- 
ginnings of XX''' century (Levyp925 1, it was not until the work 
of Shao & Nikias ( 1993 1 that they received more interest in the 
signal processing literature. The a-stable distribution is a gen- 
eralisation of the Gaussian distribution that furnishes tractable 
examples of impulsive behaviour and allows us to describe 
such behaviour by means of a small number of parameters. The 
a-stable distributions are usually defined by their characteristic 
function: 



- a=1.1 
01=1.5 



i/r(w) = exp [inw - y \w\° B,„,„} 



1 -HiySsgn(w)tan(f ) 
1 -HiySsgn(w)§log|w| 



if a 5t 1 
if a = 1 



(12) 



(13) 



Fig. 1. Probability density functions corresponding to three dif- 
ferent a-stable models. The parameters of the a-stable are 
y6 = 1,7= 1 and a = 1.1 (solid line), 1.5 (dotted line) and 1.9 
(dot-dashed Une). For an easier comparison among them, the 
p parameter of each distribution has been set so that the maxi- 
mum of all the pdfs, coincide at the arbitrary value x - 50. The 
pdf becomes more impulsive (departs more from the Gaussian 
case) as a decreases. 

where -oo </i<oo, y>0, 0<a<2 and -1 < /? < 1. The 
four parameters p, a, /3 and y uniquely and completely deter- 
mine the stable distribution. The meanings of these parameters 
are: 

1 . The parameter a is called the characteristic exponent and 
sets the degree of impulsiveness of the distribution. For 
a = 2 the distribution corresponds to the Gaussian distri- 
bution and, as a decreases the distribution gets more and 
more impulsive. Another particular case is when a = 1 
and /3 - 0, that corresponds to the Cauchy distribution. 
For a i (0, 2] the inverse Fourier transform of ijjiyv) is not 
positive-definite and hence is not a proper probability den- 
sity function. 

2. The parameter fi is called symmetry parameter and deter- 
mines the skewness of the distribution. Totally symmetric 
distributions have (3 - Q, whereas yS = ± 1 corresponds to 
totally skewed distributions. 

3. The parameter y is called scale parameter. It is a measure 
of the spread of the samples from a distribution around the 
mean. When a = 2 we get the Gaussian case and then y = 
cTq/^, where crc is the dispersion of the Gaussian. 

4. The parameter /V is called location parameter and basically 
corresponds to a shift in the x-axis of the pdf. For a sym- 
metric (yS = 0) distribution, p is the mean when 1 < a < 2 
and the median when < a < 1 . 

As it can be seen, the first case in eqs. il2\ and (I13> . a 1, has 
exactly the same expression as eq. Q. For simplicity, through 
this paper we are not going to consider the case a = 1, as it 
corresponds to a single point (of zero measure) in the interval 
(0,2]. 
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Going back to the expressions in Sect. |2 we see that eq. 
held when 1 < 77 < 3, that is, when < o- < 2. As A: 
and Q.^ are positive in eq. and < a < 2 we have that 
7 > 0. By eq. we know that p - \. Therefore, eq. Q is 
exactly the characteristic function of an a-stable distribution 
with maximum positive skew. 

The fact that a population of point sources that are dis- 
tributed in intensity following a power law and that are ob- 
served with a pencil-beam instrument produce an a-stable dis- 
tribution of deflections is very convenient. The or-stable repre- 
sentation off'ers several advantages: 

1 . Simplicity: a non-Gaussian distribution that follows eq. ( I13> 
can be completely described by only four parameters, in- 
stead of an infinite number of moments. 

2. Mathematical justification: a-stable distributions include as 
a particular case the Gaussian distribution, and share with 
it many desirable properties. First, they satisfy the gener- 
alised central limit theorem which states that the limit dis- 
tribution on infinitely many i.i.d. random variables, possi- 
bly with infinite variance distribution, is a stable distribu- 
tion (Feller [T966> . Therefore, the use of (T-stable distribu- 
tions is strongly justified from the theoretical point of view, 
as they are able to describe a wider range of data which 
might not satisfy the classical central limit theorem. In sec- 
ond place, a-stable distributions have the stability property: 
the output of a linear system in response to or-stable inputs 
is again a-stable and various aspects of linear system theory 
developed for Gaussian signals extend directly to the case 
of signals with c-stable distribution. For more information 
on the mathematical foundation of stable distributions, see 
Samorodnitsky & Taqqu ( 1994i and references therein. 

3. Ubiquity: a-stable can be shown to be the limit distribu- 
tion of natural noise processes under realistic assumptions 
pertaining to their generation mechanism and propagation 
conditions (Nikias & Shao .1995,) . They agree with empir- 
ical data extremely well in so diff'erent situations as noise 
in telephone lines, atmospheric noise, radio networks, radar 
systems, financial time series, etc. Even in cases in which 
there is not a strong theoretical or physical evidence that 
an expression such as eq. (I13> holds, a-stable representa- 
tion still provides a good modelling of many processes. For 
example, later in this work we will show that the a-stable 
model works well even when the source counts do not fol- 
low a pure power law, or when the power law is cut at a 
certain flux limits. 

There is another advantage in the a-stable formulation: they 
have been thoroughly studied in the literature, and their prop- 
erties are well understood. Until recently the a-stable distribu- 
tions were generally avoided for two main reasons: first, the 
probability distribution has not a closed form in real space 
(except for the particular cases of the Gaussian, Cauchy and 
Pearson distributions). This greatly hamper the development 
of statistical signal processing techniques such as maximum- 
likelihood and Bayesian estimates. The second one is that the 
non-Gaussian a-stable distributions have infinite variance (and, 
in some cases, as we have seen in the case of EPS and 77 > 2, 
a ju pai-ameter which is not finite), and it was considered that 



a-stable distributions cannot be physical. But the same objec- 
tion applies to a very well-known process, the white noise, but 
it is however universally used in all fields of science and en- 
gineering. The variance of the theoretical white noise is not 
finite, but this doesn't stop scientists to use it as an accurate 
and useful model for real, finite processes. In the same way, 
a-stable distributions have been shown to provide excellent fit 
to a very wide class of processes observed both in the natural 
world as well as many artificial systems. In the last few years a 
great deal of effort has been carried out in the signal processing 
field to overcome the two drawbacks above mentioned. Now 
there is a plethora of available methods to perform statistical 
inference on a-stable environments. In this work we will focus 
on the application of existent techniques for a-stable parameter 
extraction in order to obtain optimal estimators of the param- 
eters describing the differential counts of the EPS population, 
namely the slope rj and the normalisation k. 

4. Point source parameter extraction using 
a-stable distributions 

After the general discussion of a-stable distributions presented 
in the last section, let us now return to the particular case of the 
P(D) and its a-stable distribution posed by eqs. Q to (II 1> . 

According to equations (|8j to (II 1> . the usual parameters 
describing the differential counts of the EPS population are di- 
rectly related with the parameters of the a-stable distribution 
of observed deflections. In particular, using equations (|8j and 
(I10> we have that 



Tj — a + \, 



k - y- 



'"+'r(2ii)r(2i±2)sin(f ) 



(14) 



(15) 



Therefore, it suffices to estimate the parameters a and y to di- 
rectly estimate 77 and k. Over the past years a number of ef- 
ficient estimators for the parameters of a-stable distributions 
have been developed. Unfortunately, most of them consider 
only the estimate in the case of symmetric a-stable distribu- 
tions (y6 - 0), due to the fact that is the most common case in 
many signal processing applications. On the other hand, very 
recently Kuruoglu ( 12001 ) introduced a number of density pa- 
rameter estimators for skewed a-stable distributions. The sim- 
plest estimators are based on the following idea: let us consider 
an a-stable distribution with parameters a, (3, y and and de- 
note it by Sa(J3,y,fj)- Xi, i - 1,.. .,N is the sequence of 
data, it is easy to show that very simple manipulations of the 
data can be performed in order to produce centred, deskewed 
or symmetrised sequences, respectively: 



^3* + X2k+\ 

M-T 



yS,[2 + 2«]r,0 



1 + 2" 

- X^k + - 2'^"'X3J:_2 

~ 5„(0,4y,[2-2'/«]A7) 

- X2k - X2k-\ 

~ Sa{Q,2y,Q). 



(16) 



(17) 



(18) 
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In the previous equations, the symbol ~ means equahty in dis- 
tribution and therefore they must be regarded as exact equa- 
tions, not approximations. The only caveat is that they work 
only if the samples X„^, X„i_i, etc are independent random vari- 
ables. Due to the correlations introduced by the beam at small 
scales, this is not true if the samples that appear in each indi- 
vidual summation are neighbours. Fortunately, at scales larger 
than a few arcmin, where the effect of the beam is negligi- 
ble, statistical independence is satisfied. To ensure the validity 
of equations il6\ to (I18> it is enough to randomly shuffle the 
data before operating. Since in this paper we do only a first 
order statistical modelling -we deal with the ensemble of data 
rather than a time or space series- this has no other effect on 
the methods we present here apart from guaranteeing the good 
behaviour of the previous equations. 

Once the distribution is conveniently centred, deskewed or 
symmetrised, the techniques for symmetric ff-stable parame- 
ter estimate can be applied. Kuruoglu (2001 1 describes several 
groups of techniques adequate for such task: fractional lower 
order moment (FLOM) methods, logarithmic moment methods 
and extreme value methods. Other useful methods are based on 
the study of the empirical characteristic function of the data. 



Kuruoglu (1200 U studied the comparison between the different 
techniques and showed that both FLOM and logarithmic meth- 
ods are very efficient in general. In this work we are going to 
use the logarithmic method, as it is easier to implement. 

4.1. Logarithmic moments estimators 

Let X be a set of data distributed following an a-stable ~ 
Sa(J3,J,0). Let us define the logarithmic moments of the dis- 
tribution 



4.1 .1 . Logarithmic estimator for a 

Apply centro-symmetrisation as given by eq. dlSt to the ob- 
served data to obtain transformed data. Estimate L2 and then 



L2_l 



-1/2 



(25) 



4.1 .2. Logarithmic estimator for /3 

Once a has been estimated, obtain a distribution with ju = 
(for example centring as in eq. il6\ ). estimate L2 and then 



\0\ 



2 



L2 



1/2 



(26) 



Then, estimate |y6| using eq. (I24> . If centring was applied, it is 
necessary to transform the resulting yS by multiplying by (2 -H 

2")/(2 - T). 

According to eq. (|9|l for the particular case we consider in 
this work, fi - \. We provided here the expressions for the 
estimation of \Q\ and fi for the general case, but in the following 
we will make use of our knowledge of the true value of fixing 
it to its theoretical value /? = 1 instead of trying to determine it. 

4.1 .3. Logarithmic estimator for y 

Assume again that = (or make it centring the distribution 
as in the previous case), estimate L\ and hence 



y = cos(0) exp ([Li - (/tq] a -H lAo) . 



(27) 



Li =E[log|X|], 
L2=Ef(logm-E[logm])2]. 



(19) 
(20) 



where E is the usual estimator operator. It can be shown 
(Kuruoglu .200U that 



Li = (/roll - - | + - log 
\ a I a 

(I 1 \ 02 
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cos 

where i//k are the values of the polygamma function 

tpk-i = ^iogr(x) 



(21) 
(22) 



(23) 



that takes values i/zq = -0.57721566..., iff] - 7t^/6, 4^2 = 
1 .2020569 . . ., etc, and is a dummy parameter 



e = arctan^jStan^^jj. 

This leads to the following estimators: 



(24) 



Take into account that if centring was applied, y should be cor- 
rected by multiplying by 1/(2 + 2"'). 

The estimate of the location parameter is a tricky issue. 
Although there are some proposed methods, they usually show 
problems of convergence and applicability. Fortunately, the de- 
termination of jj. is not necessary for estimating the parameters 
T] and k and therefore we will obviate this issue. 

Provided with the estimators (I25> . (I26> and (I27> we are in 
position of determining the parameters 77 and k. 

4.2. Point source parameter extraction from non-ideal 
power law number counts 

In the previous sections we have showed that an extragalactic 
point source population whose number counts follow an ideal 
power law leads to an or-stable PiD) distribution. However, a 
pure power law like the one in eq. (^ is an idealisation with 
no physical meaning. On one hand, eq. ([0 diverges when S 
goes to and, on the other hand, the power law extends to in- 
finite fluxes, which is not physically realisable. From the point 
of view of astronomy, it is not possible to find galaxies of arbi- 
trarily high flux and, if we are willing to avoid Olber's paradox, 
a minimum flux has to be imposed as well (at least as long as 
ri>2). 
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In a real observation there must be a minimum and a maxi- 
mum flux Smin and Smax, respectively. This leads to a truncated 
power law 

fO if S<S,„i„ 
n(S)^l kS-'^ if S„i„ <S < S„ax (28) 

i if S > S,nax 

From the point of view of modelling point sources in a 
given sky map, if the area of the map and the number of sources 
(galaxies) are finite, the maximum flux Smax can be safely con- 
sidered as infinite, being the probability of finding an extraor- 
dinarily bright source negligible. 

An analogous situation can be found in the example of 
white noise. The theoretical white noise power spectrum is flat 
in all the range from to oo, whereas in reality it is not possi- 
ble to find such a process, but something similar with certain 
cuts that depends on several factors such as the data size, sam- 
pling, etc. Nevertheless, white noise is a very good model for 
instrumental noise and many other very well known examples. 
In a similar way, we expect that the or-stable model will be a 
good one to describe the P(D) distribution originated from a 
truncated power law like in (I28> . at least for sufficiently 'well- 
behaved' cases. 

Intuition supports the choice of a-stable distributions as a 
fair approximation to the true P{D). When observing a finite 
sample of a theoretically infinite process, if the sample is big 
enough and the process is sufficiently well-behaved, we expect 
the general model to be an adequate description of the sam- 
ple. Observing a finite number of galaxies implies a cut like in 
eq. ( I28> . that is, that we do not observe infinitely bright nor in- 
finitely faint galaxies. In any case, the basic shape of the P{D) 
distribution, an asymmetric bell-shape with a positive tail (i.e. 
the impulsive behaviour), should be preserved. The extremes 
of the tails will not be completely realised, but the basic shape 
should be kept over a certain range, orders of magnitude in size 
if the ratio Smax I Smin is great enough. 

In order to test the validity of the c-stable approximation 
as well as the performance of the logarithmic moments esti- 
mators introduced in the last section, we performed exhaustive 
numerical simulations, reproducing the observation of typical 
truncated power law-distributed point sources detected through 
a Gaussian beam in a variety of cases. We have found that the 
a-stable model is a very good approximation when the non- 
Gaussian tail of the distribution is allowed to be well-realised. 
By 'good approximation' we mean that the parameters and 
k can be estimated with errors below 5% by means of the log- 
arithmic estimators presented above. The goodness of the a- 
stable model depends on the following factors; 

1 . Sample size: there must be enough data samples to permit 
the non-Gaussian tail to appear and the logarithmic esti- 
mators to work. The simulations show that ~ 10^ or more 
samples are more than enough to work safely. Note that a 
typical 512 x 512 pixel image satisfies this condition. 

2. Flux limits: in order to clearly observe the tail of the distri- 
bution, the range of fluxes of the point sources must extend 
over an interval big enough. In other words, there must exist 




Alpha stable distribution 

Simulated Point Sources 



0.995 1 1.005 1.01 1.015 

D 

Fig. 2. Comparison of the P{D) distribution function of simu- 
lated point sources with an a-stable distribution fitting the his- 
togram data. The parameters of the EPS simulation were rj-l.l 
(a = 1.2), Smin - 10"^ (in arbitrary units), Smax = lO-' (in 
the same arbitrary units), Npix - 2048 x 2048, Ns - Npi^ and 
FWHM = 3.0 pixels (corresponding to a true value of the 
normalisation k = 3.014 x 10"** in the chosen arbitrary unit 
system). The normalised histogram of deflections is shown by 
means of a dotted line. The logarithmic moments estimators 
applied to that simulation give the estimates a = 1.184 and 
y = 2.521 X 10"^ (corresponding to an estimated value of the 
normalisation k - 3.051 x 10"'*, in the chosen arbitrary unit 
system). Using these estimates, the corresponding a-stable dis- 
tribution with jS = 1 is shown using a solid line. The position 
of the dotted line has been slightly shifted to the right in order 
to make clearer the plot. 

a relatively few galaxies much brighter than the vast major- 
ity of low-flux ones. A dynamic range S^ax/Smin ^ 10"* 
suffices to guarantee the good behaviour of the tails. This 
condition is satisfied in CMB observations as well as in 
many other observations at different wavelengths. 

3. Slope of the differential counts: as the distribution becomes 
more and more non-Gaussian, the tails grow and there are 
necessary more and more galaxies to 'map' these tails. This 
means that for a values near to 2 it is very easy to get the 
shape of the distribution with few points, whereas when a 
decreases one needs a much higher number of points and a 
wider flux dynamic range in order to realise the tails. The 
results given in the two previous points are valid for a 6 
(1, 1.5); for higher a values it is much easier to realise the 
tails and the a-stable approximation holds even for smaller 
number of data and tighter flux cuts. On the contrary, below 
a = 1 the a-stable approximation is less accurate unless the 
data size and the flux cuts grow accordingly. 

4. Beam size, pixel size and number density of sources. The 
combination of the first two factors define the so-called ef- 
fective resolution element of the experiment. The intuitive 
idea behind the effective resolution element is that, due to 
the smoothing produced by the instrumental beam and the 
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finite size of the pixel, the fluctuations in the images are 
dumped below a certain scale that is given by the beam and 
pixel sizes. If the beam is sufficiently small, the effective 
resolution element corresponds to the pixel size, but if the 
beam size is bigger than the pixel size (as is usual) the ef- 
fective resolution element is related to the coherence scale 
of the fluctuation field. For a Gaussian beam of width (Jb 
the coherence scale is roughly equal to cr^ (Rice[l954}, and 
the area cr\ defines the effective resolution element. In the 
case of a circular Gaussian beam cr/, is related the the Full 
Width Half Maximum (FWHM) by FWHM = 2o-,, V21n2. 
The a-stable approximation is shown to work optimally 
when there is approximately one source per effective res- 
olution element. When the number density of sources is 
lower, the field is undersampled and bigger numbers of 
data are needed to realise the tails of the distribution. On 
the other hand, very faint galaxies, whose number den- 
sity is much higher than one per effective resolution ele- 
ment, do not practically generate intensity fluctuations. In 
this latter case, the fluctuations are dominated by sources 
whose number density is just below the limit of 1 source per 
effective resolution element, i.e. the brightest undetected 
ones. Hence, the method is very sensitive to objects able 
to generate intensity fluctuations, that is, the source popu- 
lation which dominates the number counts at fluxes such 
that the number density is ~ 1 per effective resolution ele- 
ment or lower. Note that this limitation is not exclusive of 
the Qf-stable model, but it is shared by every other method 
that works with the P{D) distribution (Scheuer [T974t . The 
reason is that the level of ~1 source per effective resolu- 
tion element (in the sense of coherence scale explained be- 
fore) roughly determines the r.m.s amplitude of the distribu- 
tion whereas fainter sources only add some Gaussian noise 
(Franceschini et al. .1989. Barcons .1992j . Moreover, it is 
not possible to distinguish between one source per effective 
resolution element and an infinitely large number summing 
up the same total flux of an individual source. 

Therefore, for a e (1,2] the a-stable is a good approxima- 
tion, provided that the number of data and the flux limits of the 
galaxies take reasonable values that are easily satisfied in usual 
Astronomical cases. 

As an example, a 2048 x 2048 pixel simulation was per- 
formed in which point sources were distributed following a 
truncated power law with flux limits Smin - lO"-', Smax - lO-' 
(in arbitrary units) and slope 77 = o- -H 1 - 2.2, and then 'ob- 
served' with a Gaussian beam of FWHM - 3.0 pixels. The 
mean density of sources in the simulation was one per pixel. 
The resulting distribution PiD) is shown in Fig.|2| The real his- 
togram of the deflections produced by the EPS simulation (the 
solid line in the figure) is compared with an a-stable distri- 
bution with the parameters extracted from the EPS simulation 
using the logarithmic moment estimators in eqs. i25\ . i26\ and 
J27> . The agreement between the two curves is very good and 
the relative errors in the determination of the parameters a and 
7 (or, conversely, rj and k) are ~ 1 % for both parameters (see 
figure caption). 



Let us summarise the results of this section. We have pro- 
posed a set of very straightforward estimators to extract the pa- 
rameters 77 and k which characterise the differential counts of a 
generic point source population. These estimators, described in 
equations i25\ to (I27> . are based on the logarithmic moments of 
ff-stable distributions and are specifically designed to deal with 
non-Gaussian, asymmetric pdf s such as the one generated by 
point sources in the sky. From the theoretical point of view, 
these estimators are efficient (Kuruoglu, 2001) whereas 'clas- 
sical' analysis based in ordinary moments (mean, variance, 
skewness and so on) is not reliable since they do not converge. 
Moreover, the estimate of the parameters a and y, directly re- 
lated to the parameters 77 and k, is direct and computationally 
very fast, since it only needs the calculation of two moments Li 
and L2 (eqs. (ED and The analysis of a 2048 x2048 pixels 
takes around one second in a PC with a XEON 2.0 GHz proces- 
sor. An approach based on 'classical' moments would require 
the calculation of an higher number, > 3, of moments, their 
comparison with a precalculated set of values given by a cer- 
tain model and, finally, finding the best parameters by means 
of a fit (and so on). The a-stable assumption works well for the 
case of truncated EPS distributions, provided that the number 
of data is large enough and the ratio between the cuts S„,ax/S 
allows the tails of the distribution to be correctly realised. By 
using logarithmic estimators the parameters 77 and k can be es- 
timated with very small relative errors (~ 5%) for a wide range 
of 77 values. As it happens with other existent methods on the 
PiD) distribution, this method works optimally when the aver- 
age number of sources per effective resolution element is ~ 1 . 
This means that when we estimate the parameters 77 and k of the 
differential source counts we are, actually, estimating the pa- 
rameters of the source population which dominates the counts 
in the ffux interval around the S value corresponding to ~ 1 
source per effective resolution unit. 

5. Point source parameter extraction from 
non-ideal power law number counts: realistic 
galaxy population models 

In section 14.21 we have estabUshed the performance and the 
range of applicability of estimators based on the a-stable mod- 
elling of the P{D) distribution for truncated power laws, that are 
a fair approximation to the observed number counts. Provided 
that reasonable conditions of applicability of the estimators are 
satisfied, the results are valid for diverse fields of Astronomy, 
including the study of the X-ray background and the modelling 
of unresolved point sources at radio wavelengths. 

In this section we will go one step further and consider 
state-of-the-art realistic number counts models. As an example, 
we will focus on the application of the a-stable distributions to 
microwave observations and, in particular, to the images that 
the future ESA's Planck mission will produce. 

The study of real microwave images differ from the study of 
the simulations in the last section for two main reasons. First, 
number counts of extragalactic sources do not follow a pure 
power law distribution, but a more complex behaviour that de- 
pends on the emission properties of galaxies, i.e. their energy 
spectra, as well as their local densities and redshift evolution. 
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The power law distribution is only a first order approximation 
to the real one. Second, microwave images contain not only 
EPS signal, but also CMB radiation, other Galactic and ex- 
tragalactic foregrounds (synchrotron, free-free, dust emission 
and Sunyaev-Zel'dovich effect) and instrumental noise. All the 
above has to be taken into account in a realistic analysis. In this 
section we will deal with the real number count distribution of 
the sources. Some hints about how to deal with signal mixtures 
will be introduced in next section. 

5.1. Extragalactic point sources at Planck frequencies 

To test the efficiency of c-stable distributions in estimating 
the relevant parameters, 77 and k, of source number counts in 
CMB maps it is better to rely on realistic cosmological evo- 
lution models for sources. The relevant source populations 
at microwave frequencies are "flat"-spectrum compact radio 
sources, selected at cm wavelengths, and galaxies whose emis- 
sion is dominated by dust, i.e. high redshift spheroids and low 
redshift starburst and spiral galaxies, observed in the far-IR 
bands. By exploiting all the available data on extragalactic 
sources coming from surveys at cm and far-IR wavelengths, 
T98 presented a phenomenological evolution model which al- 
lowed to predict source number counts in the whole frequency 
range around the CMB intensity peak. A thorough study on 
source contributions to the intensity fluctuations of the CMB 
was also presented in that paper. As discussed in Sect.^ the ob- 
servations of the microwave sky provided by NASA's WMAP 
satellite (Bennett et al. "2003b) and the surveys coming from 
VSA and CBI experiments, strongly support the predictions on 
number counts of EPS discussed by T98, at least up to fre- 
quencies V ~ 40 GHz. Given that "flat"-spectrum compact 
sources are the dominant population in this frequency range, 
we may confidently rely on the T98 model for simulating point 
sources in the sky up to v ^ 100 H- 200 GHz. On the other 
hand, at frequencies v > 300^00 GHz many new data have 
been published since 1998 and most recent evolution models fit 
better than the T98 model the available data on source counts. 
These recent models (e.g., Granato et al. 2001. ,2004i show, 
in particular, a steeper slope of the differential counts of EPS 
at 5 ~ 10 - 100 mJy and for 300 < v < 900 GHz, where 
the contribution of high-redshift spheroids show up (see, e.g., 
PerrottaEnOD- On the other hand, at fluxes S < 0.1 mJy all 
models must converge to a sub-euclidean slope in order to not 
exceed the integrated far-IR background. Thus, in all those flux 
ranges where there are no great changes of slope, the power- 
law approximation still holds and the c-stable method can be 
effectively applied. This is the case of all (or almost) all Planck 
channels given the current flux limits for source detection fore- 
seen for the mission (see Vielva et al. 2003" where all the fore- 
grounds have been taken into account). On the other hand, the 
method could reveal not easily appliable, or not appliable at 
all, to the P{D) data that will come from the future surveys of 
the Herschel mission. In fact, current estimates on the sensi- 
tivity limits foreseen for Herschel (Negrello et al. I2004> show 
that, inside the flux range where the method could be able to 
recover the parameter k and 77 -i.e., for fluxes at the level of ~ 1 



source/beam, see Sect. 4.2-, EPS counts should show a sudden 
upturn, due to either the strongly positive k-correction of dust 
emission spectra and to the strong cosmological evolution of 
high-redshift star-forming spheroids. 

In this first application of the method we still adopt the orig- 
inal T98 model since we are mostly interested in testing the 
method rather than using it in a specific observation scenario. 

5.1 .1 . Radio sources 

Radio loud AGNs (radio galaxies, quasars and BL-Lacs) are 
expected to dominate the counts in Planck LFI channels at 
fluxes S >1-10 mJy. At frequencies around 30 GHz the typ- 
ical values for the power law slope are 77 ~ 2.0-2.15 (Taylor et 
al. l2003l Mason et al. 12003^ at fluxes 10 < 5(mJy) < 300. At 
higher fluxes, the data coming from classical radio surveys at 
cm wavelengths show typical slope greater than the Euclidean 
one (i.e., 77 > 2.5. On the other hand, at lower fluxes the 
power law index should keep ~ 2.0-2.2 down to fluxes of 
S ~ a few ju7y where it has to break down to lower values, 
for not exceeding current limits on the integrated extragalactic 
background (see, e.g., Haarsma & Partridge 1998 and refer- 
ences therein). The number of expected detections at the 30 
GHz Planck channel, based on the T98 models and using the 
Mexican Hat Wavelet detection technique (Vielva et al. 2003}, 
varies from ~ 1800 (when the emission of the rotational dust 
is taken into account in the simulations) to ~ 2700 (when it is 
not). Evidently, the number of detections depends on the flux 
detection limit attainable by the chosen technique. 

5.1 .2. Dusty galaxies 

Both 'normal', i.e. spiral-like, and active galaxies show dust 
emission that quickly dominates over the radio emission at 
wavelengths shorter than a few mm. From a theoretical point 
of view, the physical processes that govern galaxy formation 
and evolution are poorly known, but there is evidence of strong 
cosmological evolution in the far-IR/mm region, particularly 
for early type galaxies (see Granato et al. 1200 II and references 
therein). Therefore, it is not easy to model number counts of 
these source populations. SCUBA, MAMBO and IRAM sur- 
veys are rapidly providing a great amount of data in this par- 
ticular energy domain and all these data are guiding the predic- 
tions on source counts and related statistics by means of phe- 
nomenological as well as physical evolution models (Toffolatti 
et al.T998" Guiderdoni et al.T998' Granato et al."200! Rowan- 
Robinson 2001 I. Anyway, all these models predict that number 
counts of EPS are dominated by dusty galaxies at v > 300 
GHz. The number of these sources detectable by Planck is vari- 
able, depending on the emission properties of the cold dust, 
on the cosmological evolution of sources and on the capabil- 
ity of detection techniques. Current estimates by Vielva et al. 
(29Syi) based again on the T98 model with galaxies whose po- 
sitions in the sky are Poisson-distributed, predict the detection 
of ~ 12700 (85 % completeness level) point sources in the 857 
GHz Planck channel. 
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Fig. 3. Predicted integral counts in the ten original Planck chan- 
nels. The solid line shows the total number counts in the simu- 
lations (including both radio selected and far-IR galaxies). The 
dotted line is the best-fit single power law model for each chan- 
nel. The fit is performed using only number counts of bright 
sources down to the flux at which there is approximately one 
source per effective resolution element, that is, the flux interval 
in which the statistical analysis of the P(D) is sensible to the 
number counts law. 

5.1 .3. Total counts 

Taking into account the mixture of the different types of galax- 
ies and all the observational and theoretical constraints on 
them, it is evident that the number counts can not be described 
by a single power law as in eq. Q. Even a truncated power 
law as in eq. j28> is not a correct representation of the number 
counts. The real number counts need a slope 77 that changes de- 
pending on the flux range considered. In some cases, the data 
can still be approximated by the sum of two or more popula- 
tions whose number counts can be described by simple power 
laws with different slopes. 

In spite of this, the a-stable model can still be useful to ex- 
tract information about the general behaviour of the dominating 
point source population from the analysis of the total P(D). In 
order for this affirmation to be true, the shape of the number 
counts curve should be similar to a power law at least in the 
ffux range of the dominant^ point sources, and the PiD) should 
be not much different from an ff-stable . Let us see if any of 
these conditions are satisfied by Planck observations. 

Figure |3] shows the logN-logS curves for the ten Planck 
channels according to the T98 model. The solid lines show the 

^ By 'dominant' we mean the sources that produce the greatest con- 
tribution to the number counts in the flux interval that corresponds to 
number densities ~ 1 per effective resolution element, as explained in 
SectlOl 



total number counts in the simulations (including both radio se- 
lected and far-lR galaxies). The figure shows that the behaviour 
of the number counts does not correspond to a single power 
law. As discussed before, this is not surprising, given that the 
spectra and the evolution properties of sources dominating ra- 
dio and far-lR source counts are quite different. However, a 
power law fitted to the number counts (straight dotted lines), 
only at fluxes greater than the one where there is approximately 
only one source per effective resolution element, does not de- 
part much from the true curve. The slopes of these fits lead 
to best-fit slopes 77 ~ 2.1 - 2.7. The best-fit slopes and nor- 
malisations calculated for each channel are shown in table 
According to Fig. |3] The best-fit slopes could be considered 
as a good approximation to the true slope of the source popu- 
lation which dominates the counts in the relevant flux interval 
(see Sect. 14. 2> in each case: i.e., "flat"-spectrum radio sources 
in the low frequency channels, dusty galaxies in the high fre- 
quency ones. 

As mentioned before, in many cases the data can be approx- 
imated by the sum of two or more populations whose number 
counts can be described by simple power laws with different 
slopes. If one of these populations happen to dominate over 
the others in the flux range corresponding to ^ 1 sources per 
effective resolution element, the a-stable model will be a fair 
approximation able to determine the parameters corresponding 
to such population. If there is not a dominant population, but a 
number Np of roughly equally rich populations with different 
slopes, the situation becomes more complex and the a-stable 
model will not be correct in general. However, as A^,, increases 
the situation tends again to a-stability, due to the generalised 
central limit theorem (Sect. |3}- In this latter case, the a-stable 
will allow us to determine the average parameters of the mixed 
source populations. 

5.2. Validity of the a-stable model 

A very simple way to see if the a-stable model is valid for the 
Planck point source simulations is to fit the number counts to a 
single power law with certain fit parameters a^„ and kfn, then 
to apply the logarithmic estimators in eqs. j25l l, ilbi and i27\ 
to realistic Planck EPS simulations in order to obtain some es- 
timated aiog„, and A:/,,^,,,, and then check the agreement between 
the two sets of values. 

We performed realistic point source simulations using the 
T98 model. The simulated sky maps were filtered using the 
beam sizes (circular Gaussian approximation) and the pixel 
scales of the 9 Planck channels. The number counts of the T98 
model were fitted to a straight line in order to obtain the values 
of the parameters a fir and k ji, that better fit the model in the flux 
range of interest. Taking into account the results of section RT^ 
the fit was performed using the number counts at fluxes S > S i, 
where S 1 is flux over which there is approximately one source 
per effective resolution element. The parameters a and k were 
then estimated for each channel using logarithmic moments es- 
timators. The results can be seen in tabled Table^shows that 
the validity of the a-stable model is better for some channels 
than for others. The relative errors in the parameter a are only 
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Table 1. Comparison between the point source parameters a 
and k extracted using logarithmic estimators (denoted by the 
'login' subscript) and the best-fit power law (denoted by the 
'fit' subscript) for the nine frequencies of the T98 Planck EPS 
simulations. The normalisation k is expressed in units of Jy" 
pixel"^. A good agreement between the fit and the logm values 
indicates that the a-stable approximation is valid as a means 
for estimating the slope and the normalisation of the differen- 
tial number counts of the dominating point source population. 
Relative errors are calculated as r.v = 100 x \xi„g,„ - Xfi,\ /x/n, 
X representing either a or k. 
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Fig. 4. Comparison of the histogram of the distribution of de- 
flections due to extragalactic point sources in the 30 GHz 
Planck channel (solid line) and an a-stable distribution (dot- 
ted line). The a-stable distribution has the a and y parameters 
that were estimated by means of logarithmic estimators on the 
EPS data: a = 1.24 and y = 7.02 x 10-^ The EPS map is 
expressed in tsJ jT units. 



a few percent. The agreement in k is worse because equations 
M5\ and (I27> are very sensitive to the estimated value of a and 
therefore relatively small errors in a can lead to big errors in 
k. For the case of the 30 GHz channel kfn and kiogm agree up 
to a 12%, and the agreement in the 44, 70 and 100 GHz is 
even better, whereas for the worst case (217 GHz) the diff'er- 



ence between the two values rises to a 52%. For the 857 GHz 
channel the relative difference between kfi, and kiogm is 15%. 
Therefore, we can conclude that the a-stable is a good model 
for the T98 sources, and that the logarithmic moments estima- 
tors can recover quite well the most representative value of a 
for the fluxes greater than the flux over which there is approxi- 
mately one source per effective resolution element. This value 
varies depending on the channel, and corresponds in general to 
a few tens of mJy in the Planck case. The estimation of k is less 
reliable in general, but the method can be still useful, however, 
as a means to constrain its value. 

The results of Table^show that the agreement between the 
fitted a and the value estimated with the logarithmic moments 
is worse for the case of the intermediate channels. This result 
is easily explained taking into account the previous discussion 
about the physical nature of the galaxies at microwave frequen- 
cies. Whereas the low-frequency channels of Planck (30, 44 an 
70 GHz) are dominated by radio-selected point sources and the 
high-frequency ones (545 and 857) are mainly populated by 
far-lR galaxies, at the intermediate frequencies the mixing of 
two main different populations distorts the shape of the P(D) 
distribution, making the a-stable approximation less valid. As 
an example of this, in Figs.|3and|5lthe histogram of the deflec- 
tion distribution due to the EPS as two of the Planck channels 
are compared with the a-stable distributions generated using 
the corresponding estimated aiogm and kiogm- Figure |3 corre- 
sponds with the 30 GHz channel whereas Fig. [S] corresponds 
with the 217 GHz case. For this latter channel, the basic shape 
of the P{D) is very well approximated by the a-stable, but small 
differences in the region of the tail, barely visible to the eye, are 
enough to hamper the validity of the a-stable model. On the 
contrary, for the 30 GHz case in Fig.|4lthe agreement between 
the two curves is good along all the curve, and specially good in 
the tail, where the most important information is stored, lead- 
ing to a good estimation of the slope a and the normalisation k 
of the radio source counts. 

6. Point source parameter extraction in presence 
of noise and other signals 

We have shown that the tools based on a-stable analysis are 
applicable to astronomical data such as those that will be ob- 
served by the low and high frequency channels of Planck. Even 
when the mixture of two roughly equally dominant different 
populations makes impractical to determine the normalisation k 
of the mixture, it is still possible to set constraints to the global 
slope of the number counts. But the previous discussion was 
restricted only to the emission of point sources. Apart form the 
EPS population, the microwave sky contains emissions from 
many other astronomical sources ('foregrounds'), CMB radi- 
ation and instrumental noise. Therefore, it is not possible in 
general to observe the point sources independently from the 
other components. Though the addition of these other signals 
does not modify the a-stable nature of the point sources P(D), 
it affects the way we can learn the parameters a and y from 
the data. The estimate of these parameters can get very com- 
plicated in presence of all these components. In the following 
we will consider the simplified (albeit very commonly found in 
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Fig. 5. Comparison of the histogram of the distribution of de- 
flections due to extragalactic point sources in the 217 GHz 
Planck channel (solid line) and an a-stable distribution (dot- 
ted line). The a-stable distribution has the a and y parameters 
that were estimated by means of logarithmic estimators on the 
EPS data: a = 1.67 and y = 0.34 x 10 The EPS map is 
expressed in AT/T units. 

practise) problem of a mixture of a Gaussian random process 
and EPS. This is the case when instrumental noise (generally 
modelled as Gaussian white noise) is added to the EPS signal 
by the detector. In the case of CMB observations, that would 
be the case as well in 'clean' sky areas (that is, not affected by 
foreground contamination). We shall comment about this in the 
next section. 

6.1. Parameter estimation in Gaussian and a-stable 
mixtures 

Let us consider the case of a mixture of EPS 'confusion noise' 
and a Gaussian signal. This Gaussian signal can be instrumen- 
tal noise or even the CMB itself if is has not been previously 
separated from the data. The Gaussian distribution is a particu- 
lar case of the c-stable family with o- = 2. The mixture of two 
a-stable distribution with different a index is not an a-stable 
distribution. Therefore, the simple logarithmic moments esti- 
mators presented in Sect.l^are not valid. 

The characteristic function of the mixture of two inde- 
pendent random variables can be expressed as the product of 
the characteristic functions of the two original variables. That 
means that the characteristic function of a mixture of an a- 
stable distribution Sa(J3,J,f^) and a Gaussian N{0,cr) is 

^mixiw) = exp l^ifiw - y B„.,„ H- ^cr^>v^| ' (29) 

where B„a is defined in eq. ( I13> . Applying centro- 
symmetrisation as in eq. ( I18> we have 

iA^,./w) = exp [-2y I wr - o-W\ . (30) 



The parameter extraction in the case of mixtures with char- 
acteristic functions such as in eq. ( I30> is a very difficult problem 
that has not been totally solved yet. It has been studied by How 
& Hatzinakos (.1998j . who presented some basic ideas based 
on the empirical characteristic function (ecf) 

1 ^ 

«Am(w)^ -^e'-'^C"'. (31) 

The ecf is a complex random variable and its expected value 
coincides with the true characteristic function of the distribu- 
tion when the x{m) samples are i.i.d. How & Hatzinakos ( 1 998» 
described two types of methods based on the ecf to perform the 
estimate of the parameters a, y and cr in eq. ( I30> : minimum 
distance methods and moment-type methods. We tested both 
kind of methods and found that the moment-type ones suffer 
from problems of stability in the particular case under study. 
Therefore, we will focus on the minimum distance method. 

6.1 .1 . Minimum distance metliod 

In the minimum distance method, the estimate of the parame- 
ters - (a, y, cr) is obtained in the optimisation process 

|i/'Af(w) - i/'0(vv)| W(w)c/w, (32) 

CO 

where W(vv) is an appropriate weighting function. For exam- 
ple, the choice W(w) - exp(-w^) allows the integral ( I32> to be 
solved by means of Gauss-Hermite quadratures, which is com- 
putationally convenient. Thus, the estimate of the parameters 
is reduced to a minimisation over three parameters. An interest- 
ing possibility appears when one of the parameters, generally 
cr, is a priori known. Then the minimisation gets much more 
simple. 

The minimum distance method is as far as we know the 
best choice present at this moment in the literature, but it may 
present some problems. In particular, the choice of the weight- 
ing function lV(w) is arbitrary. The choice W(w) = exp(-w^) 
is useful from the computational point of view but it may mask 
important features of the ecf. We are currently working in the 
elaboration of more robust methods, that will be the scope of 
future works. For the moment, let us restrict the discussion to 
the minimum distance method. 

6.2. Some simple tests 

In order to test the minimum distance method we performed a 
set of simulations of point sources whose number counts fol- 
low a truncated power law, adding a Gaussian white noise to 
the beam-convolved sources. The size of the simulations was 
Npu - 1024 X 1024 and the cuts of the truncated power law 
were set so that S„,a_JSmm - 10^ in arbitrary units. The beam 
was set to FWHM = 3 pixels. For the sake of clarity, the simula- 
tions contains one source on average per pixel, that corresponds 
with ^ 1 source per effective resolution element (as defined in 
section l4~2t . given the value of the FWHM chosen. 

We have tested the performance of the method as a function 
of the relative contribution of each component of the mixture. 
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In order to do that, the slope of the EPS law was set to rj -2.2 
{a - 1.2). Fifty simulations were performed for each case. We 
fixed the contribution of the Gaussian noise to cr^ = 1 (in arbi- 
trary units) and let vary the value of the minimum source flux 
Sfnin- This has the effect of varying the rms of the EPS contribu- 
tion, i.e. the js of the ff-stable distribution. In other words, we 
fix the second term in the exponential of equation |30| and vary 
the first one. As described in section|3l we write the scale pa- 
rameter of the Gaussian contribution to the characteristic func- 
tion as yc - o-q/2 (where era is the Gaussian noise dispersion). 
Then we can study the performance of the method as the ratio 
7s lyc varies. Results are shown in Fig.|6l 

The method obtains its better performance when the contri- 
butions of the ff-stable distribution and the Gaussian noise are 
comparable in the characteristic function. Inside the interval 
-1 < \ogiys Ijc) ^ 1 the relative errors in the determination of 
the three parameters are of a few percent. When the Gaussian 
contribution dominates (\og{ys lyc) < -1) the estimator 'sees' 
only a Gaussian contribution and tries to adjust it to the three 
parameters by considering that the two terms inside the expo- 
nential in eq. ( I30> are identical Gaussian exponents. Therefore, 
it wrongly gives a a value of 2, underestimates the true jc by a 
factor ~ 50% and greatly overestimates js ■ 

When the point source population overdominates the 
method tends to assign the central parts of the distribution to 
an almost inexistent Gaussian and to fit the remaining tail to 
an ff-stable that is more impulsive than the real one, producing 
lower values of a than the true ones. The js parameter is well 
established in this region, but the estimated yc gets artificially 
high. 

In many real cases, the level of the Gaussian contribution 
is a priori known. For example, the noise level of an instru- 
ment is generally well known in most experiments. If that is 
the case, the minimisation ( I32> can be simplified by fixing the 
value of cr (that is, yc). Then the minimisation is performed 
in the two-dimensional parameter space (a, y). In this case, the 
estimates of the a-stable parameters improve significantly in 
the low ys /yc regime. As an example, in the previous test 
experiment for ys = 7.48 x lO"-' (corresponding to the point 
logiys /yc) = -1.82 in figure|6}, the estimates of the parame- 
ters when cr was unknown are a = 1.75±0.22,/s =0.09+0.11, 
whereas if the value of cr is fixed the parameters are much bet- 
ter estimated, a = 1.2 ± 0.1, ys = 7.54 x 10"^ + 0.91 x 10^1 
Our tests show that if cr is well known it is safe to apply the 
method until ratios \og{ys /yc) - -2.2, with few percent rel- 
ative errors in a and y. Below that threshold, the performance 
drops quickly. 

In the intermediate ys /yc regime the improvement is not 
so spectacular: as an example, for ys - 0.472 (corresponding 
to the point log('ys /yc) - -0.025 in figure|6|l, the results were 
a = 1.2 + 0.02, ys - 0.476 ± 0.02 when cr was unknown and 
a = 1.2 ± 0.01, 75 = 0.466 ± 0.006 when it was known. 

In the high ys /yc regime the knowledge of cr gives few 
information and the results are practically the same when it is 
known than when it is not. Of course, if the EPS largely dom- 
inate the mixing it may be a better idea just to forget about 
the minimum distance method and to apply directly the log- 
arithmic moment estimators of section 0] as the data can be 
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Fig. 6. Performance of the minimum distance method as a func- 
tion of the relative contribution of the mixture components. The 
top panel shows the relative error in the determination of the 
Gaussian dispersion crc- The relative error in the determination 
of the a-stable parameter y is shown in the middle panel. The 
bottom panel shows the relative error in the determination of 
the index a. Statistical error bars are shown in each case. 

regarded as an o'-stable contaminated with a small amount of 
noise. To test this idea we just applied the logarithmic mo- 
ment estimators to our simulations. The results show that if 
logiys /yc) 5 0.6 the logarithmic moments fail, as expected, 
but over that threshold they work very well. As an example, for 
ys - 7.5 (corresponding to \og(ys /yc) - 0.57) the minimum 
distance method obtained a - 1.27 ± 0.47, ys - 6.48 + 3.46, 
whereas direct application of the logarithmic moments gives 
a = 1.19 ± 0.01, ys = 7.24 + 0.12. At higher ys/yc ratios 
the validity of the logarithmic moment estimators is even more 
justified: fory^ = 118.6 (corresponding to \og(ys /yc) = 2.37) 
the minimum distance method obtained a - 0.55 + 0.80, 
ys - 22.0 ± 22.6, whereas direct application of the logarith- 
mic moments gives a - 1.19 + 0.01, ys = 109 + 4. A small 

negative bias ( 5%) seems to appear in the determination of 

ys for high ys /yc ratios, probably due to the presence of the 
Gaussian interference. 

Taking all the previously said into account, the reccomen- 
dation is to use the minimum distance method when the ex- 
pected \og(ys /yc) < 1 and just the logartihmic moments esti- 
mator for higher cf-stable / Gaussian ratios. Fixing the value of 
the Gaussian contribution with a priori information allows us 
to determine the a-stable parameters safely up to log(ys /yc) ~ 
-2, whereas if cr has to be determined the safety region finishes 
around logiys /yc) ~ -1- 

We also performed simulations varying the slope rj (a) at a 
given ys /yc (both contributions were set so that the contribu- 
tion to the characteristic function of both components are the 
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same). The estimates for the three parameters a, y and cr have 
low relative errors (~ 10%) around the values a ~ 1.0-1.6. 
Far from this region, the estimates become worse, specially in 
the low-ff regime, where bigger simulations or larger dynamic 
ranges (the ratio S „,ax/S min) are needed in order to 'map' the a- 
stable tail correctly. On the other hand, when a gets too close 
to 2 the two distributions are almost Gaussian and is very hard 
to distinguish them. 

7. Some considerations about EPS parameter 
estimation at Plancl^ frequencies 

As mentioned before, the microwave sky as observed by any 
of the many experiments that are now being carried out or in 
preparation is a mixture of CMB, galactic foregrounds, EPS, 
galaxy clusters and instrumental noise and systematics. The 
full study of the EPS' confusion noise in such a complicated 
mixture is out of the scope of this work and will be addressed 
in future studies. However, some considerations about the ex- 
pected problematics and possible results for the future Planck 
mission can be hinted in base on the results of previous sec- 
tions. 

Galactic foregrounds will be a major problem for the ap- 
plication of the techniques presented in this work. Foreground 
signature is nor Gaussian nor ff-stable, and therefore neither 
the logarithmic moment estimators nor the minimum distance 
method for Gaussian plus a-stable mixtures are expected to 
work well under strong Galactic contamination circumstances. 
Neither will do any other of the existent methods for the P{D) 
study. 

On the other hand, in areas of the sky where the Galactic 
contamination is low there would be present basically EPS, 
CMB and instrumental noise (as well as some secondary CMB 
anisotropics due to SZ effect). CMB is Gaussian or near- 
Gaussian distributed, and therefore in such sky areas the min- 
imum distance method presented in section 16.1.11 can be ap- 
plied, just considering the mixture of CMB plus noise as a sin- 
gle Gaussian component with joint effective variance cr^jj = 

•^CMB '^n- Such approach should success for the channels 
where the EPS signal is more relevant, that is, the low fre- 
quency channels (30 and 44 GHz) where radio sources are 
dominant and the high frequency channels (545 and 857 GHz) 
where dusty galaxies are dominant. 

Unfortunately, precisely at those channels the foreground 
contribution is expected to be very high. In particular, at 857 
GHz will be practically impossible to find a single patch of the 
sky where the Galactic dust emission is subdominant with re- 
spect to CMB. Foreground-free areas of the sky could be found 
outside the Galactic plane in the intermediate Planck channels, 
where EPS are expected to be so faint that they will be under 
the detection threshold of the minimum distance method. 

A more interesting possibility appears if the foreground 
emission can be somehow removed from the data before the 
analysis. As mentioned in the Introduction, there are sev- 
eral different statistical component separation methods that are 
able, given some amount of a priori knowledge about the sta- 
tistical properties of the different components, to extract them. 



These methods in general are not able to deal with unresolved 
point sources, and normally a residual map containing noise, 
EPS and some leftovers of the components that have not been 
perfectly separated is obtained as a by-product of the sepa- 
ration. Note that most component separation assume that the 
non-separable noise is Gaussian, which is not true in presence 
of EPS. This may introduce errors in the separation of the fore- 
grounds. An interesting possibility that will be studied in the 
future is to introduce our knowledge on the a-stable nature of 
EPS confusion noise in the component separation method. 

If a perfect separation was possible, the residual map would 
contain just instrumental noise and EPS. In that ideal case 
and for the T98 EPS model and the noise levels expected for 
Planck, the ratios between the point source and the Gaussian 
noise contributions to the residual map would be \og{ys Ijc) - 
2.77, 3.09, 3.31, 2.05, 0.75, 0.30, 0.17 and -0.12 for the 30, 
44, 70, 100, 143, 217, 353, 545 and 857 GHz channels, re- 
spectively. According with the results in section |6j such ratios 
should allow the methods presented in this work to obtain the 
ff-stable parameters of the P{D) with small errors for all the 
channels. 

The real case will not be so optimistic as the one just men- 
tioned, but not so pessimistic as the case in which all the fore- 
grounds are fully present. The study of the performance of the 
a-stable methods after the application of a component separa- 
tion method -such as the Maximum Entropy method (Hobson 
et al. ll999i) . for example- will be carried out in a future work. 

8. Conclusions 

In this work we introduce the formalism of a-stable distribu- 
tions as a useful tool for the modelling and the statistical study 
of number counts of undetected point sources in astronomical 
images. When the number of faint point sources in an image 
is large, the unresolved EPS contribution creates a 'confusion 
noise' that can be studied to obtain information about the par- 
ent EPS population. 

We have shown that when the differential number counts of 
point sources follow a simple power law the characteristic func- 
tion of the resultant distribution of intensity (or temperature) 
fluctuations is exactly an a-stable one. In Sect. |3l we briefly 
review the many definitions and properties concerning the a- 
stable formalism. The a-stable model allows us to model non- 
Gaussian, impulsive distributions in a flexible way and with a 
reduced number of parameters. The mathematical foundation 
of a-stable analysis is well established and useful results such 
as the generalised central limit theorem are available. The a- 
stable distributions include the Gaussian distribution as a par- 
ticular case. 

We have shown as well that, under reasonable conditions, 
the a-stable model is well suited to describe point source popu- 
lations following truncated power law number counts. We also 
have shown that even if the number distribution in flux is not a 
pure power law the a-stable model may be useful to estimate 
the parameters of the dominant source population. In particu- 
lar, the a-stable model has been proved to be useful and very 
efficient in recovering the k and 77 parameters of the EPS num- 
ber counts foreseen for the highest and the lowest frequency 
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channels of the Planck surveyor mission, where point sources 
give the dominant contribution to CMB fluctuations. As the 
EPS number counts in these channels are strongly dominated 
by flat-spectrum radio sources (at low frequencies, LFl chan- 
nels) and by dusty galaxies (at high frequencies, HFI channels), 
the method discussed here has proved helpful in determining 
the main parameters of the differential number counts of both 
populations. 

The a-stable model allows us to design efficient and com- 
putationally fast estimators to extract the physical parameters 
of the EPS populations from PiD) distribution. In particular, 
the logarithmic moments estimators are able to determine the 
slope T] and the normalisation k of the EPS population with rel- 
ative errors of a few percent over a wide range of conditions. 
The method is shown to work optimally in the case in which 
there is approximately one source per effective resolution ele- 
ment. 

Even when the EPS signal is mixed with Gaussian inter- 
ference, it is possible to estimate its parameters using the em- 
pirical characteristic function, given by eq. (13 U . We suggest a 
minimum distance method that is able to deal with the presence 
of Gaussian noise in a wide range of values of the ratio js Ijc- 

Furthermore, the method uses all the information content 
of the data, taking into account bright sources as well as very 
faint ones which contribute to the confusion noise distribution. 
Anyway, as discussed in section l4~2l the a-stable model works 
optimally in the case of approximately one source per effective 
resolution element and, thus, its maximum efficiency is reached 
at fluxes corresponding to this source density. An interesting 
possibility is to complement the information extracted with this 
technique with other methods. For example, by detecting and 
counting bright sources it is possible to obtain the slope rf and 
the normalisation k in the high flux range. These quantities can 
be compared with the ones obtained by this method for study- 
ing the differences among source populations which dominate 
the counts at high and intermediate fluxes. The way to proceed 
should be detecting first the bright resolved sources and, once 
removed these sources from the data -for example masking the 
area occupied by them-, analysing the PiD) of the unresolved 
sources by means of the statistical estimators introduced in this 
work. The statistical study of the P{D) allows us to go below 
the detection threshold. The flux limit these methods can reach 
depend strongly on the characteristic of the experiment: in ab- 
sence of noise and other foregrounds, the theoretical limit is 
the flux corresponding to approximately one source per effec- 
tive resolution element. In presence of noise and foregrounds, 
this limit will be higher and must be determined for each par- 
ticular case. 

As mentioned before, the a-stable model works well even 
there is the presence of two or more source populations with 
different slopes of the number counts, provided that the depar- 
ture from a simple power-law is not extreme in the relevant flux 
range. Then, techniques such as the logarithmic moments esti- 
mators are able to estimate a single slope that corresponds to 
the one of the dominant source population. It is possible, how- 
ever, to refine the results. A method similar to the minimum 
distance method presented in Sect. 16.1.11 can be conveniently 
modified to include more than one different a-stable compo- 



nents in order to determine the parameters k and 77 of two or, 
possibly, more source populations. 

Work to obtain an optimal technique to deal with a-stable 
mixtures is currently in progress. 

In this work we considered that the spatial distribution 
of the sources in the sky is uniform. However, the sources 
are expected to show some degree of autocorrelation. Source 
clustering will produce a broadening in the P(D) distribution 
(Barcons 1992). This effect should be taken into account in 
a future work. We did not take into account the presence of 
Sunyaev-Zel'dovich (SZ) effect in the CMB maps. The signa- 
ture of SZ clusters is expected to be similar in nature to the 
signature of EPS^ and, thus, the formalism of this work could 
be applied to it. At frequencies v < 217 GHz the SZ effect will 
have a negative contribution to the maps, producing skewed 
P{D) distributions with negative tails. Several ways of discrim- 
inating between unresolved point sources and SZ clusters in 
CMB maps have been studied by Rubino-Martm & Sunyaev 
(.2003t . 

The previous conclusions can be applied to different 
fields in Astronomy including the X-ray background, radio 
Astronomy and, in general, to all the observations in which it is 
interesting the study of the statistical properties of undetected 
point sources. 

The potentialities of the a-stable modelling go further than 
the designing of estimators such as the ones presented in this 
work. The a-stable model allows us to use techniques present in 
the signal processing literature for achieving a complete prob- 
abilistic description of the data that can be used for more ambi- 
tious goals such as Bayesian estimation, denoising, etc. Future 
works will explore these interesting possibilities. 
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